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Chapter 1 


DsTTRODUCTION 

Whenever a viscous fluid enters a stationary duct^ 
a layer of the fluid sticks to the wall of the duct creating 
a velocity gradient in the transverse direction of flow* In 
order to satisfy the mass balance, transverse velocity conponents 
towards the duct-centre and conseqiaent acceleration of the 
core surrounding the duct centreline are set up* The cross- 
flow and acceleration proceed till the fluid reaches a station 
in the axial coordinate beyond which the velocity profile 
remains constant and the cross— flow vanishes* The region in 
the duct between this station and the entrance is known as the 
hydrodynamic ent^^^ region* In this flow development region the 
momentum of the fluid increases as the velocity distribution 
becomes less uniform, and viscous dissipation of the energy 
occurs in establishing the velocity pixjfile* Hence the 
pressure gradient in this region is higher than that in the 

fully developed region* Entrance- region heat and mass transfer 

v 

coefficients are also larger than those for a developed flow, 
owing to the presence of cross flow and smellier botJndary-layer 
thickness# The calculation of these coefficients requires 
detailed knowledge of the velocity field in the entrance^region ■ 
of the duct* 

Though circular ducts are common, non-circular ducts ; 
may occur, for example, in nuclear reactors and in conpact 
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heat— exchange equipment* Since nuclear reactors use low 
Prandtl number (and hence lass viscous) liquid metals/ the 
flow development length is significant* In a compact heat- 
exchange equipment the entrance length would constitute a 
significant portion of the total duct-length* 

The hydrodynamic development of flow, in the entrance- 
region of tubes and ducts, has been a subject of study from as 
early as 19 22 when Schiller [l] presented an analysis* Even 
for laminar conditions, the flow development in the entrance 
region does not permit an exact solution for any shape of the 
duct cross— section '» The difficulty in the analysis is traced 
to the nonlinear nature of the inertia terns which appear in 
the equations of motion* In the past, various approximate 
methods of solution have been developed to provide the recjuisite 
information relating to the velocity and pressure fields* 
These methods may be brought together \inder four general 
classes x-zhich x^^ill be explained subsequently* 

In the first class are methods which linearize the 
inertia terms. With the usual boundary-layer assurtptions of 
negligible axial momentum diffusion, and negligible dependence 
of pressure on the transverse coordinates, Langhaar [2 ] 
linearized the inertia terms of the Navier— Stokes equations, 
and solved the same for a steady flow in the entrance-region | 
of a straight circular tube to get a family of velocity profiles' 
defined by Bessel functions* His results are in close agreement 
with the e:q>erimental results of Nikuradse [s]* Such linearized 
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solutions with some refinement or the other are presented for 
a circular tube by Slezkin [ 4 ], and for a paraillel-plate channel 
by Sparrow, et al * [s]* Later Lundgren, et al* [6]^ with the 
usual boundary-layer assumptions, devised a general analytical 
method for determining the entrance pressure^drop for ducts of 
-arbitrary cross sections# Their analysis has two drawbacks# 
First, it does not predict the actual entrance-length# Second, 
since the pressure-drop calculation is based on the momentum and 
mechanical energy equations, knowledge of the exact velocity 
field is essential# 

McComas [ 7 J extended the work of Lundgren, et al# [ 6 ] 
for the determination of hydrodynamic entrance-lengths from the 
fully developed velocity profile# Conparison of his results 
with the detailed analysis of development region for circular 
tube and parallel -plate channel [8,9], rectangular ducts [lO], 
and annular ducts [ll] indicates that his predictions of the 
entrance-lengths are far smaller than those predicted by [s— ll]# 
Also the centreline velocity and the axial pressure gradient are 
in substantial error# 

Flemming and Sparrow [i2 ] generalised the method given 
in [S ] for ducts of arbitrary cross section# They express the 
entrance- region velocity in terms of a deviation frcm the fxilly 
developed velocity# This deviation is determined by a least 
squares method satisfying the wall boundary conditions# They 
also use the questionable assumption that the pressure drop in 
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the entrance-^region from the momentum equation is equal to that 
from the mechanical energy equation* Though their method of 
solution is applicable to all cross-sectional ducts^ results 
have been presented only for triangular and rectangular ducts* 

In the second class are methods which divide the 
entrance- region into two parts, one near the entrance of the 
duct where a perturbation solution of the boundary-layer model 
is assumed, and the other far away from the entrance where a 
perturbation solution of the fully developed velocity profile 
is assumed* These two solutions are then patched at some 
appropriate axial location-* Schlichting [l3] determined the 
approximate velocity profile in the entrance-region of a parallel 
plate channel using this approach* Atkinsen and Goldstein [l4 ] 
extended the work of Schlichting for a circular duct* 

In the third class of analysis, initially developed by 
Schiller [l]# the duct cross-section of the entrance region is 
considered to be composed of two regions/ a boundary— layer 
developing near the wall and an inviscid fluid core* An Integral 
representation of the momentum conservation principle is enployec 
assuming a suitable velocity distribution inside the boundary- 
layer and estimating pressure from the Bernoulli equation 
applicable to the un sheared core* For circular ducts and parall? 
plate channels Schiller used parabolic boundary— layer velocity 
profiles, whereas Shapiro, et al • [is] and Siegel [ l6 ] assumed 
modified cubic and modified quart ic boundary-layer velocity 
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profiles and Bogue [iv] assumed cubic boundary-layer velocity 
profiles and power law flow- Bogue also included the radial 
momentum term in the Von-Karman^s integral method [l3^ p» 137]* 
The r^orted velocity profiles from all these integral method 
solutions vary little with the assumed boiindary— layer velocity 
profile- Tomita [l9/20] achieved similar results for power law 
flow and by the variational method- For a circular tube 
CaiTpbell and Slattery [21 ] modified Schiller's method to take 
into account viscous dissipation of energy in the boundary- 
layer- Their solution yields much better results than those 
given by Schiller'' s method especially at larger distances from 
the entrance where an increasingly larger part of the fluid is 
in boundary— layer shear and viscous dissipation of energy is 
increasingly inportant* More information about this is 
available in [22]- 

Chen [23] extended the integral momentum method for 
developing flow in' a pipe and in a parallel-plate channel at 
low Reynolds numbers- Bhatti [24] obtained a closed form 
analytical solution for the development of centreline velocity 
and pressure-drop in the entrance— region of elliptical ducts* 
Aiongwith the usual boundary-layer assunptions he also assumed 
that the core is also elliptical and concentric with the duct 
cross-section- However, for a three-dimensional flow such as 

■ r 

the flow in the entrance- region of an elliptical duct, the 
integral method requires that the velocity profile in the 
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boundary layer as well as the variation of boundary’'— layer 
thickness along the duct be assumed* For a two-dimensional 
flow such as the flow in the entrance-region of circular tube 
or a parallel -plate channel^ only the velocity profile in the 
boundar^r-lay er needs to be assumed* Bhatti also assumes that 
at the elliptical duct inlet, the streamwise growth of the 
boundary-layer is identical to that on a flat plate at zero 
incidence* However^ while the former is influenced the 
favourable pressure gradient, the latter has no pressure gradient* 
Bhatti' s results are therefore suspect* 

This integral method has been used by Han and Coopar [25] 
for triangular ducts and by Sparrow [26] for rectangular ducts* 

In the last class are solutions obtained numerically/ 
by writing various derivatives that apoear in the Navier— Stokes 
equations in the finite-difference form and solving the resulting 
set of linear or non-linear simultaneous equations* Extending 
the work of Rouleau [27],Be)doia. and Osterle [23] solved the 
entrance-region flow developm^t in a parallel-plate dhannel* 

With the usual boundary— layer assumptions the Navier— Stokes 
aquations are parabolic enabling numerical marching in the flow 
direction* They linearised the non-linear inertia terms with 
the knoTvn upstream velocity field* Wang and Longwell [29]/ 
without any boundary-layer assumptions solved the entire Navieir— 
Stokes equations using -finite-difference technique for a 
parallel -plate channel with Reynolds number 300.* They introduced 
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stream function and vorticity in the continuity and momentum 
equations to reduce the number of equations and to eliminate the 
pressure terms* They experienced some convergence problems a.t low 
Reynolds number values* 

Schmidt and Zeldin [so] proceeded exactly in the 

same way as Wang and Longwell [29] and calculated the pressure- 

defect for circular ducts and parallel— plate channels for Reynolds 

number values of 100, 500 and i0,000* The kinks in the 

developing velocity profile reported by VJang and Longwell have 

also been observed by Brandt and Gillis [3l] * Later Abarbanel, et 

al • [ 32 ] proved that the kinks belong to the exact mathematical 

are 

solution and^not merely mathematical consequence of the non- 
physical singularities, though the boundary-layer theory predludes 
such a kink* Besides Bhatti's analysis for a developing flow 
in an elliptical duct, other analyses that consider three- 
dimensional developing flows confine themselves to the triangular 
[ 33 - 37 ], rectangular [34-39] or square [33,34,40-43] ducts. 

The foregoing brief survey of the literature indicates 
that the problem of flow development in the entrance- region of 
elliptical ducts remains to be solved satisfactorily* It is 
this problem that we solve by a finite-difference method. The 
fluid is taken to be incompressible and Newtonian, and the 
usual boundary-layer approximations are used* 



Chapter 2 


ANALYSIS 


The complete equations of motion for a laminar# 
incompressible, uniform property flow in the Cartesian coordinate 
system with the velocity corrponents (u#v#w), density P# 

pressure p, and dynamic viscosity p, are [ 18]^ 

continuity : 


ax ay az 


( 2*l) 


x~momentum 


p (u ^ + V + w ~) 

ax ^ ^ay 


az' 


ap 

a 


a^ . a^^ 


I + — 2 

^ ax^ ay 


+ ^) 
az^ 

(2.2) 


y— momentum : 


^ ax ay 




az 


ay 


2 

.a V . 

+ — 2 

ax 


sh. + 

dj^ 3.2^ 

(a *3) 


z—momentum ’ 


, aw 3w . 9w \ 

P(u + W 32 ^ 


4 . 

ax^ ay^ az'" 


a z 


( 2 » 4 ) 


For the elliptical duct (Fig. 2.l) with major axis = 2a and 
minor axis = 2b, and with origin at the centre of the duct 

cross-section, the boundary conditions are 




u(x,0 < y < b, a [ l-y Vi>^] = 


0 


3U 

3z 


(x^YfO) = 0 


(x/O^z) = O 

v(x/0 ^ Y ^ to, a [l-yVto^]^'^^) « O 

v(x>0*z) = 0 

(x,yiO) •= 0 

w(x,0 X y .< b* a [l-^yVto^] = O 

w(x,y*0) = 0 


3y 

3z 


>• 


(x>Oj|Z) — 0 

u (0,y, z) = 
p (D,y,z) = pQ 


IE 

3y 


; (Xj 0 ,'z)s: O 


II (x»y,o)= 0 


( 2^5a) 


( 2»5to) 


( 2»5c) 


( 2^5cP' 


( 2*5e) 


V (0,y, z) = 0 

w (0,y, z) = 0* 

Fbr the present analysis has been assumed uniform, although 
the effect of thcs duct on the flow propagates upstream* A lucid 
disGusSion of Various boundary conditions used for u(D,y,2) for 
a circular duct is given by Shah and London [44] * We note, 
however, that the assumption of a Uniform U^ at the duct- entrance 

xS'^vBxy Boixmon^ 
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An order of magnitude analysis^ essentially that of 
Carlson [ 45 ] « may now be applied to eqns* ( 2*l) to (2*4) • 
Taking a representative length in the axial direction X/ along 
the duct-axis, to be of the order of development length L, and 
in the transverse directions y and z to be of the order of 
maximum half— widths of the duct b and a respectively, we assume 
that L » a or h* Then, if terms of 0(a/L) or 0(b/L) or lower 
compared to those of 0(i) are neglected, eqns» (2»2) to (2*4) 
reduce to 



and the dependence of p on the transverse coordinates ceases# 

Eqns» (2*l) and (2»6) constitute two equations in the 
three velocity conponents u,v and w and are hence insolvable# 

An additional equation must be specified in order to complete 
the set"# Carlson [ 45 ] has proposed and verified that for a 
square duct one possible choice is to assume that at each point 
in the duct cross-section, the transverse velocity vector is 
directed towards the duct centreline* This provides the irequired 
additional relation between v and w* Garg [46] has verified 
that this relation can be extended to rectangular ducts of 
moderate aspect ratios* Feldman, et al* [4?] have also developed 
a relation between transverse velocity corrponents for the 
developing flow in an eccentric annular duct* Following these 
successes, we also <i®velop this additional relation on the basis 
of the following :hy.pothesis * 
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In a confined flow^ the vorticity diffuses normally from 
the duct wall owing to a balance of pressure and viscous forces 
at the wall, and the fact that pressure is irrpressed normally 
upon the wall* Also in an elliptical duct, owing to the absence 
of any sharp comers, the vorticity diffusion curves emanating 
normally from the duct-wall remain normal to all equal speed 
curves* This is analogous to the orthogonality of heat flux 
lines and isotherms in a thermal problem* Such a normal 
diffusion of vorticity is possible if and only if the equal 
speed curves have the same slope at the corresponding points# 
Since the duct wall itself is an equal speed curve r^resenting 
zero speed, all the other equal ^eed curves should have the 
same shape as the duct wall and are therefore concentric 
ellipses* 

Let the equation of these equal :^eed cuivesCFig* 2*2) 
be 



with 0 < < i# The slope of the normal to the curves 

represented by eqn* (2*7) is 

Eqn* (2*8) r^resents the slope of the tangent to the vorticity 
diffusion curves at any general point (z^y)* * Upon integration, 
it gives the equation of the vorticity diffusion curves 
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y « ^ {2.9) 

here K 2 is the constant of integration*. Bqn. (2*9) indicates 
hat the shape of the vorticity diffusion cxurves is a function 
-f the aspect ratio X = a/b* For a circular duct with the aspect 
•atio unity, Eqn. (2*9) represents a family of straight line<$ 

' = K^z which are readily recognised as the radial lines* 

The above discussion indicates that the transverse 
ralocity vector is directed towards the centre of the duct along 
:he vorticity diffusion curves represented by eqn* (2*9). Based 
)n this a third relation between v and w of the form 


V 

w 


= (S) 


a^2 


Z 

2 


( 2-10) 


Ls assumed* 

Eqns* (2«l)> (2.6) and (2«10) constitute a coirplete set"** 
iirith the boundary conditions in eqns. ■( 2*5a,b, c, d). A 
rigorous solution of this set, by the finite-difference method, 
in the Cartesian coordinate SA’-stem is cumbersome because of the 
curved boundary of the duct. The treatment for curved boundaries 
by the finite-difference method are available in [43,49]. Herein 
we transform the set of eqns. (2*l), (2*6) and (2*10) alongwith 
the boundary conditions into an elliptic— cylindrical coordinate 
system wi.th velocity components (u,v^,v^) respectively* 

The constant n surfaces are confocal elliptic cylinders, 
orthogonal to the confocal hyperbolic cylinders r^resented by 

'-I* ' ■ ■ 

Also see footnote on p.ll* 
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constant IJ? surfaces (see Appendix* a)« The foci of both the 
families are same and are situated at (c^O) in the (z^y) plane# 
It is to be noted that V = represents the duct wall* 

After some leangthy algebra (see Appendix* b), the 
complete set of equations reduces to 


u 


du 

9X 


cJ 


1 


(v, 




as 


V 

an 


= ^ 1 
p 


dp 


+ + 3iu)_ (2.12) 

c^ j a dvr 

v„. (coth n - ooth T},, tanh V) 

■ S _ ' • w ^ 

(cots + coth^ n„ tan S ) 

W 


U(X,S^TI^) = 0, 

(x.,7T/2<’7) = 0, 

If (x,0,0 < ^ < 0 , 

2lB n 

ah (xrO < ,0)= 0, 


>■ 


( 2*13) 


( 2*14a) 


(x,s#h^) = 0, 


v^(x,0,0 < h < T)^) = 0, 

v^x,ji/2^r}) = 0 , 

^ (x,0 <- S < -grO) = O# 


> 


( 2 .14b ) 



15 


. 0 , 

(x/0 < « < 7T/2, o) =0, 
gv 

^ (x,0,0 < n < v^) =0, 

3V_ 

^ (x, 71/2/ = 0/ 

u(0,^i;,T?) = ■^o ^ 

p(0) - Po ' 

where v is the kinematic viscosity of the fluid and 

j (_ (S^ ^ ^ L5)/c^)is the Jacobian of the transformation 

divided by the square of the focal length c* 

It should be noted that the Jacobian of the transformation! 
vanishes at the focus and therefore eqns# ( 2*ll) to (2*13) do 
not hold at the focal point# Elimination of this singularity 
in the f inito-diff erence procedure will be ejq^lained in the next i 
chapter# 

After non—dimensionalisation with ! 

! 

vx I 

X 5 — , ! 

(nc/2) 

r = , I 

n 

z = ^ , 





( 2.14c) 




( 2.1 4d) 


1 # 
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V 


w 


(nc/2) 




V 


(nc/2) 


V 


and 


p 


;he set of eqns* (2*ll) to (2'*14) reduces to 


A (|^ + H) + DV + CW = 0 , 


3Z 


(2.15) 


ATT xi/2/,- au . „ au^ 

AU g-^ -f A ' vV 4* W 5 -^; 




aY az^- 


V = -W/B , 

u(x,y,(7) = o , 


( 2 . 16 ) 

( 2.17) 




|~(X,1/Z) 


= 0 




|y (x,o,o <■ z < a) =0 , 


1 ^ (x,o < Y < 1 , 0 ) 

V(X,Y,<7) 

v(x,o,o < z < a) 
V(x,i,z) 

1^' (X,0 < Y < 1 , 0 ) 


- 0 , 

= o , 
= 0 ^ 
= 0 , 
= 0 . 


( 2 . 18 a) 


J 


> 


(2.13b) 


J 
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and 


vhere 


and 


W(x,y^cr) = 0 , 

w(x,o < y < 1, o) = 0 , 

|y( X, o, 0 < z < a ) = 0 , 

|y^x,i,z) =0, 

u(o,y,z) = 1 , 

p(o) = 0 ^ 

A = cosh^‘(7TZ/2) ~ cos^(Try/2), 


cot(Ky/2) + coth^ n^.tan(7ty/2) 
coth(Trz/2)'-coth^ anh(7rz/2) 


g^l/2 (71/2) sinh (tiz) 

g^l/2 (7T/2) sin(7ry) 

D = , 



(71/2) 


( 2-lSc) 


( 2*18d) 



Chapter 3 

THE FINITE-DIFFERENCE METHOD 


Eqn» (2*16) is parabolic in the axial coordinate X, so 
that the classical forward marching scheme in the axial 
direction can be enployed* Since eqns^ (2*15)/ (2*16) and 
(2*17) do not involve any derivatives of the tr-ansverse 
velocity components with respect to the axial coordinate^ 
their specification at the entrance as initial conditions is 
unnecessary in the true mathematical sense* But for the 
finite-difference method that we eiploy, their specification 
at the entrance is required to start the marching procedure* 
Though theoretically the transverse velocities can be proved 
to be infinite at the entrance singular pointy we specify 
them as zero* This choice of the transverse velocities affects 
only the region near the entrance, and its effect dies out as 
we march [SO] in the flow direction-. 

A non uniform finite-difference grid (Pigs* 3*1/2) is 
imposed on the flow field- Step sizes AX, AY, and AZ are 
tcihen respectively in the X,y, and Z directions* Backx^rard 
differencing in the marching X-direction and central 
differencing in the Y and Z-directions are enployed in 
eqn* (2*l6) whereas only backward differaicing is enployed 
in eqn* (2*15)* The finite-difference forms of eqns* (2*15)/ 

( 2*16) and ( 2*17), with no-slip condition on the duct-wall 
and linearization of the non-linear inertia terms by the 
upstream flow field are 



2 


k=1 


m2 ni| m3 m4 R 

Z 


Fig. 3.2 Finite -difference grid in transformed fcompu 
tationa!) plane. 
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(Ay) (AZ) 


j/lc 


^^1+1 1 k-'^i i k^ 

+ O,. ^,V,. . , .. ,.+C.. ,.W.. . ^ , ,.+A, h dJ jJi = o^ 

(ax) 

(3*l) 


j / Tc" i+i, j , k-^^j ,yrx+u3, i , ^ 


V. . , 
j,k i,j,k 


^^i+1/ j+l/k“^i+l/ j-l/k^ 


4- W 

+ ^j,k ^i,j.k 


2(AY) 

^ ^i+1 / j / k+l""^i+l / j / k-l ^ 
2(AZ) 




(AX) 


= ^A 




(AX) 


i+ljrJ+1/k 1+1 ^ J # k 1+1/ j“*l/k 

(Ay)^ 


^ ^i+1 , j / k+l~ ^'’^i+l / j / k'^^i+i , j / k~i ^ 
(AZ)^ 


i 3*7.) 


and 


V. ^ , 

1+1/ j/k 


W . 

1+1, j/k 


B . 


( 3 * 3 ) 


j/k 


x-zliere the subscripts i/j, and k resent location in X/Y, and 
2r directions respectively and 
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AZ 


and 


AZl 

for k = 1 ( 1 ) ml 

AZ2 

for k = ml(l)m 4 +l 

ayi 

for j =1(1 )nl 

AY2 

for j = nl( 1 )n2 

AY 3 

for j = n 2 (l)n 7 * 


Eqns» (3*l) and (3*2) need to be modified at the boundary 
lines, i«e*, for j = 1 and n7, and for k = !• Moreover, eqn* 

( 3 * 2 ) needs modification at lines where the grid size changes, 
i*e*, for j = nl and n2, and for k = ml- With the inclusion 
of the boundary conditions (2*18), eqns- (3*l) and (3*2) take 
the following foriins (see Fig* 3*2) : 

(a) along the boundary FR (j = 1, k l) : 


V, 


-1/9 i+l,2,k ^,2 

(AYl) 


^ ^i+1 , 1 , k-M "’^i-fl , 1 , k ^ 


(AZ) 


^ S,k^i+l,l,k‘^l,k 


^^i+l,l,k''^i,l,k^ 


(Ax) 


0. 


(3*4) 


and 


1,3c i/l,h 


^ ^i+l , 1 , k+l''^i+l , 1 , k-1 ^ 


2(Az) 


\,3c^ i,l/3c 


^ ^i+lr 1/ 3c~^i, 1 , 3c^ 

(Ax) 
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1+1 


(AX) 


(AYl) 


i+l^l^k+l ^ 1+1 /l/k"*" i+ltfl/k—l^ 
(AZ)^ 


( 3»5) 


(b) along the boundary PO (k = 1, j /: l) : 


^1/2 

J/1 


^ "^1+1/ j +1/ l”^l+l / j # 1 ^ 
(AY) 


W 


.1/2 




4 - I- r I V ^ ]3 V 

j / 1 ^ ^ 23 ^ ) j / 1 i+1 # j f 1 




(AX) 


= 0, 

( 3 *6 ) 


and 


J/1 1/j/l 


^ ^1+1 , j +1 / l”'^i+l/ j-1, 1 ^ 
2(AY) 


+ A . ^U. . , 




(AX) 


(Pi^^-Pi) ^^1+1, j+ia“^^i+l/ j/l'^^i+l/j-'l/l^ 

jwiAl _j ^ I > . . »i<«p ni , ^ *' ' ■■" 'f'' '' ■"■■ ■■■"'■' 


(Ax) 


(ay) 


2 


^ 2~"'^l+l,j ,l ^ 

(AZl)^ 


(3*7) 


(c) along the boundary OS (j = n?)^ 
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V . ( w -w ) 

l+l/n3,k ^^2 ^ i+l,n7,k+l i+l,n7/k 

n7,k (Aya) n7,k ^^ 2 ) 


^n7,k'^i+l,n7,k'*^n7,k 


n7/k"'^i,n7/k^ 
( AX) 


- 0 , (3-B) 


and 


I)^C w . „ , 

n7/k i,n7,k 


^^i+l,n7/k+l~'^i+l,n7,k~l ^ 


2 (Az) 


'^n7,k^iyn7/k 


( U . „ , -U . „ , ) 

^ i+l#n7/k i,n7/k 

(Ax) 


^ ^ ^^i+l,n8,k"’^^i+l,n7/k^ 


"n7,k 


(Ax) 


(AY3) 


i+3./n7/k+l i+l.^n7^k i+l.^n7^k*-l 

(AZ)^ 


( 3*9 ) 


Let us now derive the form of eqns* (3*l) and (3*2) 
applicable at the focal point j=i^k = l* 

We know that the focal point of the ellipse is a singular 
point of transformation since the Jacobiein of the transformation 
vanishes at that point* To avoid use of this point in the 
computations, a point is considered very near to the right 
of focal point F such that FP^ = 6 (Fig* 3*3)* The reason for 
not considering F^ on the line FG will be explained later* Hence 
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-t 

for the f inite^diff erence procedure^ F would be the first 
point (i,i,i). The discretised form of the axial momentum 
eqn • (2»16) at (i#l/l) takes the form^ 


■^l,l^i, 1/1 


(ax) 


W 2 ^~^^itl/l/l'*''^^i+l,l, jr^i+l/l/S^ 

+ ^11 W 

1/1 1/1/1 2(AZi) 


= -A 


1/1 


(Ax) 


^ ^^i+1 / 1 , l"" ^^^i+1 / 1 / 2"**^ ^^i+1 , 1 , 3~^^i+l / 1 / 4 ^ 
3(AZi)^ 


^ ^^i+l/2/l”^^i+l/l,l^ 
(AYi)^ 


( 3-10) 


when and are expressed by the second order accurate 

qZ 

/\TT 

forward differences [Sl]/ as the boundary condition (x,Y,o) = 0 
does not hold at F^ and hence use of central differences is not 
possible* 

Had vre considered F^ on the line FG (Pig* 3*3)# the 
discretised form of C2-16) at the new F^ would have been 


■^l/l^i/l/l 


^^i-fl/l/l‘'^i/l/l^ 


(Ax) 


^ 2 V 

1/1 i/ 1/1 


^ '“^^itl, 1/ I'^^^i+l / 2/ l’'^i+l / 3/ 1 ^ 


2(AYl ) 
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= ~A. 


1 # 1 


( Ax) 


(AZl)^ 


x-1-1 , 1/1 i+1 / 2/ 1 X4-1 / 3/ 1 i-f 1 / 4/ 1 

+ # 

3(AYl)^ 

qtt 

since the boundary condition ^ (X^O,z) =0 is not applicable 
at the new F^« It is clear from the last term of the above 
equation that the band width of the coefficient matrix (to 
be explained later) of the resulting set of linear Simultaneous 
equations would have been doubled in comparison to that using 
eqn» (3*10)* The consequence of this increase is additional 
computer time and storage requirements* 

The discretised continuity equation at is 

1/2 ^i+1 '’I 1/2 ^i+1 12 ^ i+l/l/l i/1/1 

(AYl) (AZl) (ax) 


(AZl) 


( 3 . 11 ) 


It is to be noted that no correction for AZi has been applied 
either to eqn. (3*10) or to eqn>, (3-11) and is tahen 

to be zero as 6 is negligibly small* 

Along lin^es where grid size changes take place, eqn* (3*2) 
takes the following form [see Ref* 50 , Appendix d]* 

(a) Along the line k = ml (all j except j = nl and n2) ; 


A^/ ^ V 

j/ml i,j,ml 


f II ™ U . ) 

i+l/J+l/ml i+l, 3 -l#ml 

2 (AY) 
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^a 1/ 2 


1 +lf j/ni3 1 1 + 1 / 3 /ml 1 ±+l,j/m2 


(l+e^)(AZ2) 


+ A u. . . 

J / ml 1/ j / ml 


^^i+l#j/ml "* ^i/j/ml^ 
(Ax) 


= -A. 






+ 


i+l/j+l/ml ^ 1+1/ j/ ml ^1+1/ j~l/ml^ 

(AyT^ 


26-1 ~ (i+e.) u. . , . + 

4. — 1 — -+-"*'l-/- 3/’^3 X J-+l/3/ml 1 1+1/ 3/ m2 


(1+6^)(AZ2) 


where 0, 


AZ2 

AZl 


< 1- 


(b) Along the line j = nl (all k except k = ml) ; 


< 3*12) 


A1/2 y 

ni/k l,nl,k 


(e^ 


2 1+1 / r4/ k 


+ (i--ej)u 


- u. 


2 i+l/nl/k l+l/n3/k 


) 


(l+e2),(AYl) 


+ ^ W 

^ nl/k ^l,nl,k 


^^l+l,nl/k+l “ ^i+l/nl/k-l^ 


2(AZ) 


^nl/k ^l/Hl/k 


^^i+l/Hl/k ''/^l/Hl/k^ 
(AX) 


^®2^®2^i+l/n4,k''^^‘^2^^1+l/nl/k'^^l+l/n3/k^ 

.A — ±±i — i- 4 . 

ni/k 


(AX) 


(l+e^XAYl)^ 


*‘^l+l,ni/k+i ~ ^^l+l/nl,k ^ ^i+l/nl,k-l^ 

(AZ)^ 


(3 •13) 
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where ©2 ~ 
(c) Along 

aV2 

n2/k 


AYl 

— < 1 * 
AY2 


the line j 
(U 


V 


±j n 2, Tc 


n2 (all k except k = ml) : 

i+1^ n6 , k" ^ 3 ^^1+1 , n 2, k*^ ^i+l , n5 , k ^ 

(l+e3)(AY3) 


4- W 

n2/k ^i,n2yk 


^^i+l/n2,k+l ~ ^i+l,n2^k~l^ 
2(AZ) 


“^nZ/k ^i,n2/k 


^^i+l,n2ik ** ^i,n2/k^ 
(AX) 


(AX) 


^®3^^i+lyn6,k'’^^''^3^^i+l/n2/k'^3^i+i,n5,k^ 
(H6 3)(AY3)^ 


^ ^i+ 1 , n'2/ k+l'’ ^^i +1 , n 2 , k'^^i+ 1 , n 2 r ^"1 ^ 

(AZ)^ 


(3.14) 


AY3 

where 0 ^ = < 1 • 

ZY2 ~ 

(d) At the point j = nl and k = ml : 

^®^i+l#n4,ml"^^^'^2^^i+l,nl,ml''^i+l,n3,ini^ 

4 y. 

ni^ml i,nl/mi 


(1+0 2^^^^!^ 
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nl,mi 


(l+ej)(AZ2) 


'*’ ■^nl,ml^i+l^nl,ml 


(AX) 


nl/inl (Ax) 


20 2 t^®2^i+l,n4/ml "■ ^^'*^2^^i+l,nl,mi ’^i+l/n3,ml^ 


YY-O ^ ^^-4 -L.-f rx-1 m-l y-xH mol 

1 i+l/Hli-mS 1 1 i+l/nl/mz-* 


( 3*15) 


( e ) At the po int j = n 2 and k = ml 


A^/ 2 y 

n2,ml ijrn2/ml 


^^i+1 , n6 ^ ml'' ^ 3 ^ ^ 1+1 , n 2, ml''® 3^1+1 , n5 , ml^ 

(l-+03)(Ay3) 


'^n 4 #ml^i,n 2 ,ml 


^^i+1 , n, 2 ^ m3*’^ ^ '’® 1 ^ ^i+1 , n 2, mi''® 1® i+1 , n 2, m 2^ 


(l+e^)(AZ2) 


4 - A U 

n 2, ml i, n 2/ ml 


^®^i+l#n 2^, nd.‘'®^i,n 2, ml^ 
(AX) 


= -A 


'n2/ml 


(AX) 
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(1463) (AY 3 ) 


26 ^ nO m-a “* ml r> O 

1 1+1/112# in3 1 i+lin2#inl 1 l+l/n2/rn2 

4. ^ 

(n^^) (AZ2) 

(3*16) 

Eqns* ( 3 * 1 ) "to (3*16) constitute a set of ( 3xm4xn8 4-m4— n3+l ) 
linear simultaneous equations in that many unlcnowns, viz* 


u . 

14-1 # J * k 

(for 

j = 1 ( 1 ) 

n3+l. 

k = l(l)m4) 

V 

i+l# 

( for 

j == 2 ( 1 ) 

n8. 

k = l(l)m4) 

W , ^ . ■, 

1+1, j,k 

( for 

j = 1 ( 1 ) 

n8+l. 

k = 2(l)m4) 

w. 

1+1, 1,1 

f 





and 


i'+l 


This rather large set of equations can be reduced 
to a set of only (m4(n8+l)+l) equations in that many unknowns ^ 
viz* 


and 


Ua ,• V ^ 3 ' = 1(1) n 3 tl ^ k = 1(1) m 4 ), 

' 'l+X^ ^ 


^i+1 ^ 
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by the use of a^cial momentum equation along with the integral 
form of continuity equation# Once U' s and P are knovn at 
the (i+l) location, V's and w' s can be found from the continuity 
eqn# (3«l) and the relation (3«3)# 

The continuity eqn# (2*15) integrated over the duct 
cross~section by the application of the trapezoidal rule is 


[ 




+ E 

k=2 2 


+ 


m4 “^l/k^i+l^l/k®! n3 l^i+1^ j # 1 


+ E 

k=m3 2 


+ E 

j=2 2 


A U. (l+§«) ^ A. U. . § 

nl,l H-l,nl,l^-'- ®2 n5 j,l ^2 

+ S 

4 j =n4 2 


^n 2 /l^i+l/n 2 , ^2'*'^4^ n8 ^j/l^i+l#j/l ^4 


+ E 

j=n6 2 


^n7,l^i+l^n7,l ^4 n3 m2 


+ S E A. , Uj . , 

4 5=2 k=2 ^ 

n 3 , mi^i+i , j , ml ^ n 3 m4 


+ E 
j=2 


2 


+ E E A ,U. - . 

5=2 k=m3 1 + 1 / 3 /k l 


m2 *nj,lc'^i+l/nl,k‘^'^*2’ *'nl.mi^i+l,nl,mi'l+®l+*2+*3> 

-f- E ■ ' " •* " ■ ' " * ' " ; + ' ' " ' ' ' ' •' ' 

k=? 2.,' '4 
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m4 

+ E 

]c=in3 


nS 

j=n4 


m2 
+ £ 
lc=2 


m4 
+ E 
k=an3 


n8 

+ E 
j=n6 


m2 
+ E 
k=2 


m4 
+ E 

]c=in3 


A._ ,U, 




A. 


;( §0+^0^ 


ns 

m2 

E 

E A 

Js=n4 

k=2:. ■ 

nS 

m4 

E 

E 

j =n4 

k=in3 




f. 


^n2, k^i+i,n2jrk^ §2"*'^4^ '^n2#ml^i+l#ri2^nil^ ^2'*‘^3'**^4'^^5^ 


'^n2j,k^i+l,n2,]c^^3‘^^5^ n8 m2 

2 j=n6 ]c=:2 i-+UU^ 4 


j # ml^i+ijr j , ml^ ^ 4+^5 ^ n8 m4 

+ S B v^i+1 -5 V 

2 j=n6 ]c=an3 5 


'“n? , k^ i+1 #n 7 < k ^4 ^nl , ml^i+x , n7 # ml ^ ^4‘*'^5 ^ 

: -.; — — — — — - + — 


^n7#k^i+i,n7,k ^5 


] (AZi)(AYi> 


1 a 

I S U,A*dZ dY 
o o 


slnh . (2T)^) 

5r~ 


( 3.17 ) 


^1 = 


AZ2 

AZl 


f ^-5 ^ 


Ay.2 

AYl 


where 



3 


$3 = ^ i *§2 

and 

% = 


9 


«4 


AY3 

AYl 



Chapter 4 


COMPUTATIONAL DETAILS 

For the solution of finite— difference equations given 
in Chapter 3# the following steps are followed 

1 * ^i+I j Tc ^i+l from egn« (3*2) 

and eqn* (3*17) solving a set of [m4(n8+l)+l] 
simultaneous equations (solution procedure to be 
explained later) - 

2* V. , are calculated from the continuity eqn*(3«l) 

1+1# J/K 

in a step-wise manner, starting at j = n8# k = m4 
and moving in the negative Z-direction upto k = l, 
utilising the above U. . , and the relation (S-'S)- 

iT*l # J # 

Simultaneously W. . , are calcxilated using the 

# J / 

relation (3*3)« 

3* Step 2 is r^eated in the negative Y-direction to 
j = 1* 

4* n7 k calc\iLated starting from k = m4 to 

3 W / \ 

k = 1, utilising the boundary condition (X,i,z; = O 
and antisymmetric nature of the function Bj (in the 
relation 3*3 )• This completes the solution at the 
location ( i+1 ) • 

5* Steps i to 4 are repeated till the flow is found 
to be fully developed* Mathematically this would 
occur at X -*• «>« Practically, ctxnputation is 
terminated when the centre-line velocity attains a 
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value equal to 99 % of the fully developed 
centre-line velocity which is 2#0 in dimensionless 
t erms • 

4 #1 Solution Procedure 

The system of (n+l) simultaneous equations resulting 
from eqns* (3*17) and (3*2) can be written as 


U Q 

[S] { } = { } (4*1) 

P R 

where n = m4(rTe+l), 


E O 

S = [ ] is the coefficient matrix of order (n+l) 

G H 

E = (E^,E 2 , E^) is the row vector arising 

from the coefficients of axial velocity 
coiTponent U at location (i+l) in eqn* (3*1?), 

H = is the column vector arising 

from the pressure at location ( i+l ) in the 
momenfum eqn* (3*2), 

[g] is an un symmetric banded matrix of order n 
and band-width ( 2xm4+l) arising from 
coefficients of axiail velocity coirponent U 
at location (i+l) in eqn* (3*2), 

P is the unknown pressure at location (i+l), , 

Q is the known right hand side of eqn* (3*l7), 
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and 

R = { known right hand side 

column vector arising from the terms at location 
(i) in eqn* (3»2)* 

Solution of the system (4»l) by the inversion method 
is highly expensive even on the presenb-day computers since 
the number of steps required is of order (n+l)^ and the storage 
requirement isof order (n+l) * A futile attempt was made to 
use the sparse- solver routines MA23AD, MA^BD and MA 28 CD of the 
Harwell library in order to utilise the sparsity of the matrix 
[s]* The counterparts of these routines in the NAG library 
also did not work due to some untraceable problem in the 
particular NAG library routine on our system* The algorithm 
suggested by Gupta and Tanji [52]J also did not work due to 
the first row of [S] consisting of the vector (E) being too 
full* Gupta and Tanji' s algorithm required that no row should 
be more than half filled in single precision mode and inore than 
one third filled in double precision mode* 

An algorithm was therefore developed to solve the 
system ( 4 *1) based on the partitioning of [S] into the row 
vector (e), column vector {H} , and the un symmetric banded 
matrix [g]* The banded matrix [g] is converted to an upper 
triangular matrix by Gauss elimination [53] with column 
pivoting* The same operations are simultaneously carried out 
for £H} and {R} also* Then the operated matrix takes the form 
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E 0 


U 


Q 


<a fr p 


(4*2) 


where [G'*'’] is the upper triangular matrix with diagonal 
elQnents unity and {H*} and { R*} are the operated vectors {H} 
and { R} respectively# If xve define 


and 


and 


and 


then 


D„ = , 

n n ' 


n-i 


s= H* — j] G£ (for i = n'f'l/n— 2/ » • 2/ 1 ) 

j _ j * j - j 


"n - ^ 


n-i 


_ R* - ®i, i+j ^i+j ' ^ ~ . **,2/1 ) 


and 


U, 


f = C_ — P #0 f 
n n n ' 


(4#3) 


Ui = Ci - P.D^ 


(for i = n-l,n~2/-» -^ 2/l). 


The first equation of the system (4»l) is 


n 


E E.U. = Q# 
i=l 

Substitution of eq^n# (4-3) Into eqn- (4-4) leads to 


(4.4) 


n 


( E E.C.) - Q 


n 

E 

1=1 


(4-5) 


E 
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and back substitution in the system (4*2) gives the axial 
velocity vector {U} i 

4*2 Gorrection to the uniform inlet velocity u^ 

With the assumption of unifoarm inlet velocity 

integral form of the continuity eqn* (3<»17) is 
not fully satisfied since the left hand side of (3#17) contains 
errors 0(AZ) and 0(AY) injected through the application of 
the trapezoidal rule# Hence a correction to the uniform 
inlet velocity applied such that eqn# {3#i7) gets 

satisfied at X = 0 itself# This corrected inlet velocity 
U . V i® given in Table 4#1 for various aspect ratios X# 

Of j » JC 

In factf it was observed that the pressure P increased for 
about 7 initial marching steps (i = l#2f#»#f7) when this 
correction was not applied# 

4*3 Problem with the transverse velocity c omponents near the 
focus 

Very high wiggles were observed in the transverse 
velocity compononts V and W in the irrmediate vicinity of the 
focal point of the ellipse# These wiggles were found to grow 
rather than darrp out in the marching X— direotionf This was 
found to be more pronounced at higher a^ect ratios and at 
finer grid sizes AZ and AY due to very finely spaced points 
near the focus It is clearly due to the singularity of the 
transformation at the focus# To avoid these 'wiggles^ the 
values of W near line OP (Fig# 4#l) were estrcpolated from 
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those away from the focus such that W(X,Y,0) = O, After some 
numerical experimentation^ it was found that a fifth order 
polynomial (exact fit) provided results accurate to fourth 
decimal place for all aspect ratios* 

It should be noted that this extrapolation has been 
carried out only over 25% of the real flow domain* 

4 •4 Selection of a value for 6 

A numerical experiment was conducted for the aspect 

ratio X = 4 with various values of 6 = (Fig* 3*3) viz* 

lO"'^ f 10**^ and 10~^^* Variations in the values of 

pressure and velocity coitponents were observed respectively in 

the fifth and sixth places of decimal between 6 = 10 and 

10 * For lower 6 values the variations were observed only 

in the tenth place of decimal* For the results presented here^ 

— 3 

6 was therefore taken to be 10 for all aspect ratios* 

4*5 Step sizes 

For computational efficiency but without sacrificing 
the accuracy the stqp size AX was changed systeanatically 
along the X— direction for various aspect ratios (Table 4*2)* 

It is to be noted tbat the dimensionless axaal distance is 
defined as , 

„ 9 9 

X r: ' ■ ■ ■ where = a^-b'^* 

inc/2) Uq 

For higher aspect ratios such as X = 2 and 4 the denominator 
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2 

(.TCc/2) Uq becomes large* Hence AX must be inversely 
proportional to the square of the focal length c* 

Computer codes (Appendix C) were developed for 

a) uniform AY and AZ values> and 

b) non-uniform AY and AZ values (specifically, 
three st<^ sizes in the Y-direction and two in the 
Z-direction ) • 

However since we did not have enough time to experiment with 
non— unifoimv grid code, the results presented here were 
computed using the unifojrm grid code with AY and AZ values 
as given in Table 4*1* Confutations were carried out in 
double precision mode (l5 significant digits) on the DEC-1090 
system at ITT Kanpur* 



Chapter 5 


RESULTS AInTD DISCUSSION 

Results for the cDomputad ajtlal velocity U(X/Y, z), centre- 
line velocity U(X, l,o), pressure P(x) and the pressure defect 
K(x) are presented for various values of the aspect ratio* 

In order to limit the number of plots to a reasonable number, 
the axial velocity U(X,Y, z) is shown only along the major and 
minor axes* Comparison is made with the already available 
results for elliptical ducts [6,24] and in the limiting case 
for a circular duct [so] . 

5 -1 Axial Velocity U(X,Y^Z) 

As the transformation from the Cartesian to the elliptic- 
cylinder system does not hold for X = 1 # an elliptical duct of 
aspect ratio l» 00 i was considered so as to check the computer 
programne against well known results for the developing flow 
in a circular pipe [ 50 ] « For this a^ect ratio the focus is 
very near the origin (q/b being only 0»045)* Fig# 5 »i(a) shows 
the developing axial velocity at the major or minor axis for 
X = l*00l (solid lines) and the same for a circular pipe by 
Hombeck's method [ 50 ] (dashed lines)# The^ velocity profiles 
at the major and minor axes exactly coincide which itself is 
a check of the present an'Slysis, as the velocity profile for a 
circular pipe is the same at any radial line# The small 
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differences could foe attrifouted to the aspect ratio foeing not 
exactly one. 

Having gained this confidence, higher values of the 
aspect ratio were tried* Pigs* 5-l(fo) to 5*l(j) show the 
developing flow axial velocity profiles at the major and minor 
axes for aspect ratio X = 1*1, 1*2, 1 •4, 2*0 and 4*0* Both the 
ellipt ic~cylinder coordinates and the Cartesian coordinates 
are indicated* Numerical values are listed In Tafoles to 

5 *1 ( j ) • The small kinks ofoserved in the velocity profiles along 
the major axis in the region surrounding the focus, specially 
at hioher aspect ratios, are due to the fact that the elliptic- 
cylinder coordinate system clusters more grid points near the 
focus, specially at higher aspect ratios* Prom the Cartesian 
coordinates sho^'m, it is clear that a uniform AY in the elliptic- 
cylinder coordinate system corresponds to a highly non-uniform 
Az in the Cartesian coordinate system* These kinks disappear 
when U along the major axis is plotted against the Cartesian z, 
as is evident from Pigs* 5*l(h) to 5*l(r)* 

These figures show that the difference in the velocity 
profiles foetween the present analysis and Bhatti's analysis [ 50 ] 
increases, in the viscous boundary-layer region, with the 
aspect ratio at the small X values* By assunption Bhatti's 
velocity profile in the boundary-layer r^ ion is almost linear 
at the small X. values for all ?\.** Bhatti also assumed -a value 
for the boundary-layer thickness in the developing region, and 
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it turns out to be smaller than that found in the present 
analysis • The above diffejrences increase with A. but decrease 
as the flow develops. The later is expected since Bhatti's 
velocity profile in the developing region is based on the fully 
developed velocity profile* 

5*2 Centre-line Velocity 

The centre-line velocity development for elliptical 
ducts of higher aspect ratios is shown in Pig* 5 *2(3) (solid 
lines) along^vi'ch Bhatti's results (dashed lines). As mentioned 
above the bound ary- layer thickness in Bhatti's analysis is 
small. Thus only a smaller portion of the fluid is in boundary- 
layer shear according to Bhatti'^s analysis. Consequently the 
centre-rline velocity is undej>-predicted* The difference 
between these two analyses is therefore maximum at the smaller 
X values for all aspect ratios* A® x increases# this difference 
decreases fear all aspect ratios. Infact# for A = 4 # Bhatti's 
analysis underpredicts the centre-line velocity at small X lout 
overpredicts it at larger X* 

Pig* 5.2(b) shovjrs the centre-line velocity development 
for nearly circular duct with x = l'*00i for the pros ent analysis, 
Bhatti's analysis [ 24 ] and for the circular pipe fso]. The 
present results coirpara-- better with those using Kombeek'^ s 
analysis [50] than Bhatti^'S results* While the ma:;d.mum 
difference between'the pfresent and circular pipe results [so] 
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is only 4 •2^ tnat between Bhattl'^s and circular pipe results 
[ 50 ] is as rmch as lO It should be noted that Bhatti's 

results for circular pipe do not compare well with the 
experimental data, especially at the small X values [24, Fig. 3]» 
His velocity underpredicts the experimental data just as it 
undexpredicts the theoretical results* 

5*3 Entrance Length 

The hydrodynamic entrance lengths for various aspect 
ratios are given in Table 5 .2 in terms of X and X'*^ (X’^ being 
X»(3Tc/2b)^) • Computations were terminated when the centre- 
line velocity reaches 99^ of the fully developed centre-line 
velocity, viiich is 2»0 in dimensionless terms. In the 
computations 'b' is always leapt unity while * b.* is varied to 
achieve different values for the aspect ratio* As per 
diseassion in Sec. 5*2, Bhatti^s analysis overpredicts the 
entrance-length for small X values and underpredicts the same 
for high X values. 

5 *4 Axial Pressure Drop and Pressure Defect 

The axial pressure drop P(x) and the pressure defect 
K(x) are shown in figs. 5.3(a) to 5.3(e) for various values of 
X • Values of P are also listed in Tables 5*l(a) to 5.l( j). 

The pressure' defect is the .difference between the actual 
pressure drop and the ideal pressure drep resulting from the 
a -^.cmirm-hion of fullv developed flow at X = 0 . This ideal 
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pressure drop can be easily shown to be 
dP/dX = nH-K^^TC^) 

since the fully developed velocity profile is 
U =2(1 - ( 2 /a) 2 - (y/b)2) . 

The difference between the present analysis and Bhatti's 
analysis is high in the initial development region for all 
aspect ratios# As e3g)ected from the above discussion, this 
diffcir Jice decreases as the flow develops# However as X 
increases, it vanishes for moderate aspect ratios (X = I*!, 1#2 
and l#4) but reduces to a non-zero value for higher aspect 
ratios (X = 2*^ and 4#0)# 

W^ien the flow becomes fully developed as X -*■ c=>, K attains 
a fully developed value, which for various aspect ratios is 

given in Table 5*2# Bhatti's analysis predicts to be 7/12 
while the analysis of Lundgren, et al • [e] predicts it to be 
2/3* However, it is known that the latter value is over 
predicted* The present analysis shows that it depends weakly 
on the aspect ratio* 

Prom the above comparison, it appears that Bhatti^s 
analysis does not provide good results in the initial developing 
region due perhaps to somewhat arbitrary assunptions in his 
analysis regarding the growth of boundairy layer region and the 
velocity profile in this region. However, a firm judgem^t 
must await cortparison with e:}^erimental data for developing flow 
in an oll,lptigal'' duct*.'. 
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Chapter 6 
CONCLUSIONS 

^he problem of flow development in an elliptical duct has 
been solved by a finite-difference method following transformation 
of the governing eguations of motion into an elliptic— cylinder 
coordinate syston# In order to solve the simultaneous finite- 
difference aquations which result from the above method, an 
algorithm was developed to partition the coefficient matrix and 
thus utilise the banded structure of the partitioned matrix* This 
saves considerable computer time* The comouted solution is 
compared with Bhatti's solution [24]* The differences between 
the present and Bhattl's results are more in the developing region 
n ear the entrance plane, specially within the bound arY~l^ysr, but 
they vanish as the flow develops* This is due to Bhatti's 
assumption regarding the growth of boundary-layer region and the 
velocity profile in this region that is based on the fully developed 
velocity profile*. The present results compare better with those 
for the limiting case of developing flow in a circular pipe than 
Bhatti's results* However, comparison with experimented, data, 
non-existant for developing flow in an elliptical duct, is the 
only real test* 

It should be mentioned that the transformation from the 



8 3 ^ 

■the focus of the ellipse and large near the wall :» A numerical 
method on the o'ther hand requires a non^^uniform mesh of just 
the opposite nature due to large gradient near the wall* Thus 
a highly non^-uniform mesh in the corrputational domain 
( elliptic-'cylinder syst em) should be used- The present analysis 
and comnutcr code (in Appendix G) is based on a non-uniform 
mesh taut due to lack of time results could be obtained only 
for a uniform mesh in the cora'putat ional domain* 
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■Appendix A 

Elliptic— cylinder coordinate system can be obtained 
by' talcing an orthogonal family of confocal ellipses and 
hyperbolas dLn a plane and translating them in the third 
direction, normal to the plane (Pig* A.i) * The relation 
between the ell iptio- cylinder coordinate (x,4,'T)) and the 
Cartesian coordinate (x,y, z) may be writtai as 


X = X , 

y = c sinh 71 sin f , - (A«i ) 

z = c cosh n cos f • 

The suirfaces ij = constant are the elliptic cylinders 

^2 2 

“ M + M = 1, (A,2) 

(c cosh v) (c sinh 7?)^ 


while the surfaces = constant are the hyperbolic cylinders 


(c 



cos ^ ) 


1 



(c sin 



1 * 


(A# 3) 


Equation of the duct wall is 



If equation (A, 2 ) were to equation, (A*4) for: ^ ss 



8S 



Fig.AI Elliptic- cylinder coordinate system. 
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then 


and 


a =s c cosh V,, t 
w 

b = c sinh , 
c = (a ^— , 


T)^ = tanh~^(b/a)* 


(A_*5) 



•%>pendix B 


Elliptic cylinders represented by constant tj surfaces are 

, / 

■ 2 ^ ■ 

“ ^ ~ + 5 = 1. ' (B.i) 

(c cosh v) (c sinh rj) 

Slope of the tangent to this surface at any general point (W^v) 
in the plane is given by 


dy 

^ = 


-cot ^ tanh T) • 


(b*2) 


Hence the slope of the normal to this suirface at any general 
point in the plane is given by 


dy 

^ = tan f coth 7). (B.a) 

We shall define a variable (Fig. B.i) as 

tan 0 = tan f coth T) • . (B.4) 

The relationships between the transverse velocity corrponents 
v^ and v^ in the elliptic-cylinder coordinate system and the 
similar velocity coupon ents v and w in the Cartesian coordinate 
system are (Fig. B.i) 


■it/2 
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Fig. B.1 Velocity components in cartesian and elliptic- 
cylinder coordinate system. 


FrxDm eqns« (B*5) and (B«6) it can also be seen that< 


V = cos e + sin 0, 


T) 


and 


w = ~v^ sin 0 + cos 0« 


We have the operators 


1 - as , ^ . 3 . 

dv dv 82 an 87 ' 

and 

2 as i- 4. siZ 

8 ^ 8 ^ 82 33 ; 37 * 

Prom Cb* 9) and (B*io) 


and 


since 


3- = ris a- as . 3-i/r9Z as - ai asi , 
87 ^8Sf8n 8n 8W^' ‘•an 8* a^f an-J 



and 


- H = c cosh n sin f , 


87 

8 l 


^ = c sinh n cos f • 

8 n 


(B*4) 


an 


= -tan f cosech^ f)/sec^ , 


(B,7) 

(B,8 ) 

(B ,9 ) 

(B.io) 

(B.ii) 

( 6 . 12 ) 

(B.13) 

(6.14) 

(b.ib) 


From eqn. 
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^ = sec f coth T)/ sec*^ 0. (B#i6) 

When the operators (B.ii) and (B,i2) operate on eqns. (B-7) 
and (B*3) respectively, we get 




30 , 


If = ® + ’'■n ® ali + ® 


0V“ 

“* v^yj sin 0 r~) + (■■ ^ sin 0 

^ dv 3f af ° 


90 

+ cos e ^ cos 0 


-V* sine If) 


(B*i7) 


and 


3 W r /• 

fz = 0 - Sin 0 


3h *" 37?. 


sin 0 


- Vj, cos e ||) 1^ + cos e 

90 

- sin 0 ^ 0 


3v 


V 


-V^ oose2f)||]/[(ff)2^(||)2] 


37) 


(3*18) 


On addition and sinplif ication (B*17) and (B*18) give 
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^ ^ 9w r ( j ^ 

ay 3 ^ ~ L*'cCv^ slnh 77 cosh 77 


+ sin ^ cos ^)/(cosh^ 77 -cos 2 f)l/2 


-c( cosh^ 7) -cos^ ^)V2 

dfl_ 


dv 


+ / [-c^( cosh^ 77 -.cos^ ^f)] 


(B*19) 


Hence the continuity equation (2.l) transforms to 


c(cosh^77 -cos^' |H + 1^ [v,^ cosh^ 77-COS^ 


+ (cosh^ 77 -,cos^ 5 l)^^ ]= 0 


in the (.Xt'a/rj) system* 

Successive operations of (B*ii) on u give 


(B.20) 


and 


an 

ay 


3 ^u 


ay 


ris M ajii /D 

^a« dv aTiaf-' ' * 


r ai a._ (M) ^ m ^ (a.ii)i /d 
^ af a 77 ^ay 377 an; ^ay^^ ' 


(B*21) 


(B*22) 


whea-e 


D 


r 15 - az a^ 1 
^ dv dw , 31 btt-* 


2 2 
•C CcDSh'^lf*- ^ 











Successive operaticais of (B#i 2 ) on u give 


- r ^ ^ -] /D 


and 


^ _ 

9 z '- 37? 3T? aw a w 


^2 
a u 


(B*23) 


3Z- + 111 * 'E. 24 ) 


aw aw 'az' 

Addition of eqns. (b. 22) and (B,24) with further simplification 
gives 

- ri3i ^ / r ..2/ v2« ' 2 

172 TT^ “ '■ 73 73 J / L c ( cosh n -cos^ W ) 1 

3y az aw a rr 

(B#25) 

From eqns. (B.7), (B.s), (B.2a) and (B. 23) 


'^ly =■ e + e) (f| - |2 ^)/d. 


and 


(B#26 ) 


w 


^ = (v^ cos © ~ v^ sin e) (- ^ ^ 2 ii)/D, 

dz V w an at) aw aw ' 


(B.27) 

of eqns. (B.26) and (B.27)/ with the substitutions 


and 


sin e = sin W coshn/Ccos^ W sinh^ rj+sin^^/ cosh^ t?)# 


cos e = cos W sinh n/(cbs^ W slnh^ 77 +sln^ W cosh^ 7 //, 


and further slmplificati^s ireduces to 




Substitlitlon of eqns* (B # 25 ) and (B*23) In the axial momentum 
eqn • ( 2 *6 ) gives 

fv M + V 


m 


^ c( CO sh^ r) -cos 2 )j; ) 1/ 2 


^ ^ 1 ^ 

o a3c 


a dir 


p dx , 2*1 2 „.\ * 

^ c ( cosh V •"COS f ) 


(B#29) 


Substitution of eqns* (B*?),, (B*8)*- (A*i) and tA« 5 ) into 
eqn* (2*10) gives 


v^ cos 0 + Vff sin 0 


-Vij, sin 0 + Vi^ cos 0 


= coth '7),. tan f tanh H 
w 


or 


V [co t » + coth^ 7)^ tan S ] 

''•■*** SS "»*• ..MHwi ii ' iii r «ll■l> wta ^^■ !n ■ IIII M IIII— 111 m m 

[coth '7?''-*<xath totnh 


w 


(B., 30 ) 


Boundary conditions t 

(a) The condition 9u(x»O^2)^0y = 0 


Now 1^ = “d 


f] sin f -t sinh I) cos '$ |~y t~*=^(cosh^~cos 
9fl OX' 


(B=.3i) 


fhen y = 0/, either 3! or ■n' is zero (Fig* B»i.)* 


or y) = 0« eon*. (B*3l) 
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and 

for ^ = O eqn • (B»31) gives 

(x,o, o < n < - 0* 

(b) The condition 3 u(x,y^O)/9 2 = 0 : 

Now If = c [sinh cos ^ 1^ -cosh T) sin W ||]/[c2(cosh2 

When z = 0, W = TC/2 and eqn* (B*32) gives 

H (x,V2#^) = Q. 

(c) The condition v(x,0,z) = 0 : 

This gives 

v^ (x/0 < U; < 71/2/0) = 0 , 
since 0 = 7r/2 along TT = 0 (Pig* B.i), 


and 


(X/0,0 < T? < = 0 


w 


since 0=0 along 'S =: 0 (Fig* B.i) . 

■ f 

(d) The condition 3v(x/y/0)/3z = 0 : 


av r/ \ . 

Now ^ = L V sinh T) cosh n sin f. cos f) “* 3$ ' 

, 

+ (cos^^+cosh^T?) 


I 


I-' ¥^4 ‘ * 


: . .--/a./-, .^ 4 ;, . 


n-cos^ ) ]• 
(B* 32 ) 


(B.33) 

(B* 34 ) 




95 


- v^cos ^ sinh ll)/(cosh^ 77~cos^ U?) 

f 2 2 

+ ( sinh^ 'V cos^ ^ r~ 

af| 

9v 

+ cosh^ V sin^^ ■— ~)]/[c( cosh^ tt-cos^ ^)V2 j ^ 

(B«35) 


When- z = 0, $ = Tr/2 and eqn# (B#33) gives 

I 

3v^ 

• (B.36) 

( e) The condition w(x#y#o) = o : 


This gives 

v^ {x,n/2,Ti) = 0, 

since e =71/2 along f = 7T/2 (Fig* B#l) • 
Substitution of this in eqn* (Bf36) results 


3v 


as 


f (x,7r/2,f)) = O* 


(f) The condition 3w(x/0,z)/3y = 0 ; 

Now ^ = [(sin^ S + sinh^ 77)cosh V cos S(v^ cosh V sin S 
AV ' T) 


+ v^ sirih T7 cos f )/( cosh.^lf —Oos ^ f )■ 


+ sinh Ti 






- 1 ; . 



- (cosh^ t) sin^ 


$ 


¥r) 


3v* 

sinh^ t) dos^ 1 / [cCGOsh^'n- cos^ ]. 


3^^ 


(B-37) 

For 11 s: o, eqn* (B»37) gives 


( xifO < liB < Ti/ 2,0) = cot f . (B .38 ) 

For 1= 0> eqn* (B»37) gives 


£j 

i’St. 


( X»0^0 < t) « 


= V^ GOth T) • 


(B*39) 


Substitution of eqns* (B*33^ (B‘«34) rospeGtiveiy in eqns* 

(B*38) and (B*3^ ) iresults 


3V«, 

(x,0 < t< »/2»0) = 0 
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